function [g gt gtt] = h2metric(params, t, u, udt, udtdt)
  % input -- T, ut, ux, uy at a single time step, of size GRID_SIZExGRID_SIZE

  tOut = params.metricEvaluationPoints;

  if (length(tOut) == 0)
    error('no time points supplied')
  end

  tord = params.timeOrder;  % order of time derivatives 




  % load precomputed function.

  if params.metricLinearFnInterpolation
    global nh2s7 nk2s7 na2t2 na2t4 na2t5;
    global interpgrid;

    load('pde/nflin.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5', 'z');
    interpgrid = z;

  else
    global nh2s7 nk2s7 na2t2 na2t4 na2t5;
    if params.metricMark
      load('pde/nfmark.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');
    else
      load('pde/nf.mat',     'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');
    end
  end


  % assumptions: all of the precomputed functions are on a chebyshev grid, from (0, 2). 

  global g0dd g0UU g0det;
  [g0dd g0UU g0det] = h2basemetric(params);




  % compute indices 
  tOut = sort(tOut);
  tOutI = round((tOut-params.Tbegin+params.Dt)./params.Dt); % index of requested time element
  tPass = (tOutI - tord/2 <= 0) | (tOutI + tord/2  > params.Tsteps);
  [tOut; tOutI];

  if (sum(tPass) > 0)
    tWrongI = tOut(find(tPass))
    error(sprintf('Desired time points are outside of available time range (%f, %f)', t(1+tord/2), t(end-tord/2)));
  end



  % compute which g need to be computed so that all time derivatives can be evaluated. basically add tord/2 padding around each tOut
  baseI = [];
  for i=1:length(tOutI)
    baseI = [baseI tOutI(i)-tord/2:tOutI(i)+tord/2];
  end
  baseI = unique(baseI); % index in the tOut index space

  % figure out which indices of baseI correspond to tOut that we actually want to compute
  baseOutMatch = zeros(1, length(baseI)); % index in the baseOut index space
  for i=1:length(tOutI)
    baseOutMatch = baseOutMatch | (baseI == tOutI(i));
  end
  baseOutI = find(baseOutMatch);


  % display result:
  [tOut; tOutI; baseOutI];


  % **********************
  %  compute g
  % **********************

  disp(sprintf('computing metric at order %d', params.metricOrder));

  % allocate memory for intermediate g computation:
  g = zeros(length(baseI), params.Nx, params.Ny, params.Nz, 4, 4);
  

  for i=1:length(baseI)
    ti = baseI(i);
    tic;
    g(i, :, :, :, :, :) = h2metriceval(u(ti,:), udt(ti,:), udtdt(ti,:), params);
    gtime = toc;
    disp(sprintf('computed metric at t=%f in %.1fsec', t(ti), gtime));
  end


  % **********************
  %  compute time derivatives of g
  % **********************

  gt  = zeros(length(baseOutI), params.Nx, params.Ny, params.Nz, 4, 4);
  gtt = zeros(length(baseOutI), params.Nx, params.Ny, params.Nz, 4, 4);

  d1k = d1ker(tord);
  d2k = d2ker(tord);

  for i=1:length(baseOutI)
    spliceI = [baseOutI(i)-tord/2:baseOutI(i)+tord/2];
    %gt(i, :, :, :, :, :) = reshape(conv2(g(spliceI, :), d1k, 'valid'), 1, params.Nx, params.Ny, params.Nz, 4, 4);
    gt(i, :)  = conv2(g(spliceI, :), d1k, 'valid')./params.Dt;
    gtt(i, :) = conv2(g(spliceI, :), d2k, 'valid')./params.Dt./params.Dt;
  end

  % ******************************
  %  splice out intermediate g
  % ******************************

  baseOutDel = find(baseOutMatch == 0);
  g(baseOutDel, :, :, :, :, :) = [];

end

